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Abstract 



We establish a connection between Optimal Transport Theory (see |Vi] 
for instance) and classical Convection Theory for geophysical flows [Pej . Our 
starting point is the model designed few years ago by Angenent, Haker and 
Tannenbaum |AHT] to solve some Optimal Transport problems. This model 
can be seen as a generalization of the Darcy-Boussinesq equations, which is 
a degenerate version of the Navier-Stokes-Boussinesq (NSB) equations. 
In a unified framework, we relate different variants of the NSB equations (in 
particular what we call the generalized Hydrostatic-Boussinesq equations) 
to various models involving Optimal Transport (and the related Monge- 
Ampere equation |Br t ICaj ). This includes the 2D semi-geostrophic equations 
[Hot ICNPl IBBt ICGPt ILoj and some fully non-hnear versions of the so-called 
high-field limit of the Vlasov-Poisson system [NPSj and of the Keller-Segel 
for Chemotaxis jKgllJLl lUMPg] . 

Mathematically speaking, we establish some existence theorems for local 
smooth, global smooth or global weak solutions of the different models. We 
also justify that the inertia terms can be rigorously neglected under ap- 
propriate scaling assumptions in the Generalized Navier-Stokes-Boussinesq 
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equations. 

Finally, we show how a "stringy" generalization of the AHT model can be 
related to the magnetic relaxation model studied by Arnold and Moffatt to 
obtain stationary solutions of the Euler equations with prescribed topology 
(see immiEl lScl lVMll lNi]). 

1 The Angenent-Haker-Tannenbaum model for Opti- 
mal Transport problems 

In this section, we consider the model introduced by Angenent, Haker and 
Tannenbaum |AHTj . This model was designed in order to seek the solu- 
tions of some optimal transport problems as equilibrium states of a suitable 
dynamical system that could be efficiently solved on a computer. The con- 
crete applications have been computer vision, image registration and image 
warping. 



1.1 Optimal transport and rearrangements 

Let us briefly recall some typical results in Optimal Transport Theory, such 
as the polar factorization of maps. More precisely, let D be the closure of a 
bounded connected open set in i?'', with a boundary of zero d-dimensional 
Lebesgue measure. Up to a rescaling, we assume the Lebesgue measure of D 
to be 1. Given an L? map y : D ^ W^, we call image measure of the Lebesgue 
measure on D hj y the unique nonnegative (Borel) measure fi defined by: 



f{x)fi{dx) = / f{y{a))da, (1) 

Rd- Jd 

for all compactly supported continuous functions / on i?'^. We have 
n{dx) = 1, / \x\'^fi(dx) = / \y(a)\'^da, 

Rd JRd JD 

which means that belongs to the set Prob2{R'^) of all (Borel) probability 
measures /i on such that / \x\'^fi{dx) < oo. In this space, we say that a 
sequence /i„ converges tightly to /i in Prob2{R'^) , if: 

f{x)fln{dx) / f{x)l2{dx) 

Rd JRd 



for all continuous function / on such that 

l/(^) 
x<^R'i 1 + 



sup , \[\^ < +00. 
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Given two maps y and z from D to R^, we say that they are rearrange- 
ment of each other if they define the same image measure. When ?/ is a 
rearrangement of the identity map, we say, in short, that y is Lebesgue mea- 
sure preserving. 

Next, we define the class of maps with convex potential: 

Definition 1.1 We say that an map from D to belongs to the class C 
of maps with a convex potential if there is a lower continuous convex function 
p : R —>■] — oo, +00] such that, for Lebesgue almost every point x E D, p is 
differentiable at x and its gradient Vp(x) coincides with y{x). 

Then, we get from [Brj : 

Theorem 1.2 (Rearrangements with convex potentials) 
For each L map y : D ^ there is a unique rearrangement map with a 
convex potential y* G C. The mapy* depends ony only through the associated 
measure /i defined by (CP- 

In addition, the nonlinear operator y E L"^ ^ y* E L"^ is continuous as well 
as the induced operator fi G Proh2{R^) y* G L^, with respect to the tight 
convergence. 

We get more precise results if y is a non degenerate map, in the sense 
that the pre-image of every Lebesgue negligible set is also negligible: 

Theorem 1.3 (Polar factorization of maps IBrf ) 

Let y be a non degenerate map from D to R. Then, there is a unique 
"polar factorization" y = Y o X where Y belongs to C and X is a Lebesgue 
measure preserving map of D. In this decomposition, Y is the unique rear- 
rangement y* of y in C and X is the unique measure preserving map of D 
that minimizes fj^ \X{a) — ?/(a)p da. In addition, X can be written: 



where ^ is a convex Lipschitz function defined on . 

For (much) more results on optimal transport, we refer to Villani's text- 
book [Vi]. The expression "optimal transport" comes from the fact that y*, 
among all possible rearrangements y of y^, is the unique minimizer of the 
"transportation cost" . 




(2) 
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where | ■ | denotes the Euchdean norm. The name "map with convex poten- 
tial" is due to CaffareUi [Caj . The concept of polar factorization has been 
extended to Riemannian manifolds by McCann [Mcj . Examples of concrete 
applications of optimal transport techniques to natural and computer sci- 
ences can be found in |FMMS| [HZTX] . 



1.2 The AHT model 

The AHT model is an attempt to get the unique rearrangement y* of ?/°, 
with convex potential, as the equilibrium state at t = +00 of the following 
set of evolution equations: 

% + {v ■ V)y = 0, (3) 

Kv + Vp = y, V.v = 0, (4) 

where y = y(t,x) G R'^, v = v{t,x) G W^, p = p(t,x) G R depend on t > 
and X G -D, and is a "dissipative" operator to be chosen, for instance 
K = I or K = —A. In these "AHT" equations, we denote the inner product 
in i?'' by ■ and we use notations: 

d d d"^ 

'-^^i j=l,d '^^3 j=l,d "^^j 

The boundary conditions for the AHT system ( P|^ are: 

i) the initial value of y at t = 0, y{0,x) = y^{x), 

ii) V is parallel to the boundary dD if K = I and v = along the boundary 
ifK = -A. 

Notice that neither p nor v need initial conditions. As a matter of fact, as 
K = I, the second AHT equation (jlj) just corresponds to the "Helmholz 
decomposition" of ?/ as a sum of a gradient field and a divergence-free field 
parallel to the boundary dD. The field p can be recovered by solving the 
Poisson problem: 

Ap = V.y, 

inside D with inhomogeneous Neumann condition Vp ■ n = y ■ n along the 
boundary, where n denotes the outward normal. Then, we get: v = y — Vp. 
So, we can write: v = Py, where P is a linear singular integral operator 
bounded in all space for 1 < p < +00, provided that the domain D 
is smooth enough. In the case K = —A, in a similar way we can write 



4 



Convection and optimal transport 



V = PaU, where Pa is a linear singular integral operator bounded from to 
the Sobolev space ly^'^, for all 1 < p < +00, D being assumed to be smooth. 
Thus we can write the AHT system fl3TO in a more abstract form: 

dty + {PkV ■ V)y = 0, (5) 

with Pk = P if K = I and Pk = Pa a K = -A. 



1.3 Expected long time behaviour of the AHT model 

Let us now explain why the AHT model is expected to solve the Optimal 
Transport (or rearrangement) problem, at least for a large class of data. 
First, we observe that equation ([3]) expresses, at least formally, that, at each 
time t, y{t, ■) is a rearrangement of y^. Indeed, for any smooth compactly 
supported function /, we get: 

^ f f{y{t,x)) dx = - [ {Vf){y) ■ {v ■ V)y dx 

dt JD JD 



V ■ V[/(y)] dx, 

(using the chain rule) which is zero, since v is divergence free and parallel 
to dD and is therefore orthogonal to any gradient field. (Notice that 
this calculation can be made rigorous provided that v has enough regularity. 
According to Ambrosio's recent improvement of the DiPerna-Lions theory on 
ODEs [Ml EL], it is enough that v belongs to L]^^{R+,BV{D,R^)).) 
Next, we get the following balance law for the AHT 

— j -\y{t,x) — x\^ dx = — I {v ■ Kv){t,x) dx, (6) 

which implies, at least formally, that y and v respectively belong to the func- 
tional spaces L'=^{R+,L'^{D,R'^)) and L'^{R+,Hk{D)). Here Hk{D) denotes 
the Hilbert space of all divergence free fields w{x) G R'^ for which w-Kw dx 
is finite with suitable boundary conditions {w parallel to dD as K = I or 
w = on dD as K = — A). The formal proof of ([6]) is as follows: 

4 / hy - x]"^ dx = - [ (y-x) ■ {{v- V)y) dx 

dt JD 2 JD 
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(using the first AHT equation (JSj)) 

= — [ V ■ Vi-\y — xp) dx — I V ■ {y — x) dx 

JD 2 JD 



D 



V ■ {y — x) dx 



(since v is divergence free and parallel to dD and is therefore orthogonal 
to any gradient field) 

= — [ V ■ {Kv + V(p — dx 

JD 2 

(using the second AHT equation (jlj)) 

V ■ Kv dx. 



D 

At this point, we can describe the expected long time behaviour of the AHT 
system through the following heuristics. Since Kv-v is space-time integrable, 
we first argue that, as t +00, v presumably tends to zero. Then we 
expect y to have a definite strong limit y°° in L^, which is, then, necessarily 
a rearrangement of y^. Passing to the limit in the second AHT equation 
(jlj), we conclude that ?/°° must be a gradient. Therefore, at the end of the 
process, y^ has been rearranged as a map ?/°° with a potential. Observe that 
this potential needs not being convex. This is obvious in the special case 
when is itself a map with a potential which is not convex. Indeed, then 

y{t,x) = y^{x), f(t,x)=0, 

is a trivial stationary solution to the AHT equations fl3lH|) and we get ?/°° = y^ 
as a map with a non convex potential. So, we need further assumptions on y^ 
to be convinced that ?/°° has a chance to have a convex potential. A natural 
assumption is that y^ is smooth with a positive jacobian determinant valued 
in some interval [r, 1/r] with < r < +00. Indeed, for such an initial 
condition, the AHT equations have a global smooth solution y (at least in 
the case when K = —A, as discussed in the next subsection), with a jacobian 
determinant that must stay in the same interval [r, 1/r], since y is always a 
rearrangement of y^. So if the convergence to y°°, as t — *■ +cxo, is strong 
enough, we expect y°° to have a convex potential and, therefore, coincide 
with the unique rearrangement of y^ with convex potential y* provided by 
Theorem 11.21 The results obtained in [AHTj are only partial and leave as an 
open question this issue. 
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1.4 Wellposedness of the AHT equations 

From the PDE viewpoint, it is crucial to check that the AHT system fl3ll^ 
is wellposed, which is done by Angenent, Haker and Tannenbaum in [AHT] . 
for a class of dissipative operator K including K = I. Let us briefly discuss 
the wellposedness issue in the cases K = I and K = —A. 
For K = I, the AHT system is similar to the inviscid Burgers equation, 

dty + {y ■ V)y = 0, (7) 

since Pk behaves like a pseudo-differential operator of order zero. Thus, the 
local in time existence of smooth solutions for smooth initial conditions can 
be obtained from rather standard energy estimates. It is a challenging and 
interesting open question whether the Lipschitz norm, in space, of such so- 
lutions may blow up in finite time (as it would be the case of the inviscid 
Burgers equation). In sharp contrast, in the case K = — A, smooth solutions 
are clearly global in time. Indeed, from ([3]), we immediately get that \y{t, x)\ 
is uniformly bounded by the sup norm of y^ that we denote by M° and sup- 
pose, here, to be finite. Thus, because of (jl]), according to standard elliptic 
regularity theory, the L^(PV^'^) norm of v{t^x) is controled by M° for all 
finite p. Thus, the same is true for the sup norm of Vf (t,a;). Differentiating 
(IHl) in X, we deduce that the sup norm of Vyit^x) in x cannot grow, in sup 
norm, faster than exponentially in t as soon as V?/° has a finite sup norm. 
So, there is no possible blow up of the Lipschitz norm of both v and y and, 
therefore, by a standard argument, smooth solutions must be global in time. 
(Notice that the dissipative operator K = (— A)^/^, with appropriate bound- 
ary condition, would be borderline to get such a Lipschitz estimate. In that 
case V would be a priori only Log-Lipschitz, just like the Yudovich solutions 
of the 2D Euler equations [MP].) For more details, we refer to the paper by 
Angenent, Haker and Tannenbaum |AHT] . where different kinds of operator 
K are considered. 

1.5 Interpretation of the AHT system in terms of Con- 
vection Theory 

From a Fluid Mechanics viewpoint, the AHT equations look very similar to 
the Boussinesq equations for convective flows, in particular to their Darcy- 
Boussinesq version. A classical model for Convection Theory is provided 
by the Navier-Stokes Boussinesq (NSB) equations that we are now going to 
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review with more details. Using the Boussinesq approximation, the Navier- 
Stokes equations for an inhomogeneous incompressible fluid subject to gravity 
along the Xd direction read: 

Po{dtv+{v-V)v)-uAv + \/p = y, V-v = 0, (8) 

dty + {v ■ V)y = 0. (9) 

Here, v = v(t,x) G R'^ is the velocity field, p = p(t,x) G R the pressure 
field, po > the average density of the fiuid, u the (constant) viscosity of the 
fiuid, while y has only one component in the "vertical" direction Xd, which 
is —g 6{t,x), where g is the gravity constant and 6{t,x) is the difference be- 
tween the density of the fiuid at (t, x) and the averaged density po of the fiuid. 
(Usually in Convection Theory, a diffusion term is added to the advection 
equation for y |ID] .) We recall that the Boussinesq approximation amounts 
to consider a variable density incompressible fluid for which the density vari- 
ations are sufficiently small to be neglected in all terms except the gravity 
force. This approximation is widely used for ocean and atmosphere modelling 
|Pe] . (To the best of our knowledge the justification of this approximation 
is still an open problem in mathematical Fluid Dynamics, mostly because 
of our rather poor knowledge of the Navier-Stokes equations for inhomo- 
geneous fiows, see discussions in [Lil IMaj for instance.) By neglecting the 
inertia term, or equivalently by setting po = in the NSB (Navier-Stokes- 
Boussinesq) equations, we get the simpler Stokes-Boussinesq (SB) (related 
to large-Prandtl- number Convection Theory as in [DORll^^, for instance). 
If, in addition, the diffusion term —uAv is replaced by a friction term such as 
V, we get the Darcy-Boussinesq (DB) model. We immediately see that both 
the SB and the DB equations are just particular cases of the AHT model 
f PBl) . for which the vector valued function y has only one component along 
the Xd axis. Indeed, the DB and SB models then respectively correspond to 
the choice K = I and = — A in the second AHT equation (jlj). According 
to the discussion made in subsection 11.31 we expect, for the AHT model, the 
y{t,x) to converge, as t ^ +oo, to a map y°°{x) with, hopefully, a convex 
potential. In the particular case of the convective DB and SB models, y{t, x) 
has only one component in the Xd direction, namely —g 6{t,x). Interest- 
ingly enough, the convergence toward a map with convex potential, exactly 
means, for the DB and SB models, that the density field tends to a density 
"profile" p°°{xd), depending only on the vertical coordinate Xd, and monoton- 
ically decreasing. This clearly corresponds, in terms of Convection Theory, 
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to a "stable hydrostatic equilibrium". Notice that a similar discussion can 
be found in Moffatt's paper [Moj (section 2) as a prelude to his Magnetic 
Relaxation model that we will consider at the end of the present paper. 

2 Generalized Navier-Stokes-Boussinesq equations 

The interpretation of the AHT model in terms of Convection Theory suggests 
the following "GNSB" generalization of the NSB (Navier-Stokes-Boussinesq) 
equations: 

e{dtv + {v ■ V)v) + Kv + Vp = F{x, y), V ■ v = 0, (10) 

dty + {v-V)y = G{x,y), (11) 

where y = y(t,x) G R"^ is a vector-valued function (m > 1, in practice 
m = d or m = 2d for the models discussed below), F and G are given 
smooth functions with bounded derivatives up to second order, respectively 
defined on R"^ and R"^ x R'^, e > is a scaling factor introduced to single 
out the inertia term, and is a linear dissipative operator. Depending on 
the applications in view, only the following cases will be considered: K = 
(no dissipation), Kv = v (linear friction), Kv = —Av (viscosity). 



2.1 Existence theory for the GNSB equations 

For simplicity, we consider in this subsection the domain D to be the unit 
periodic cube T'^ = R'^/Z'^, in order to avoid technicalities due to spatial 
boundary conditions. For the three possible choices of the dissipative oper- 
ator K 

Kv = 0, Kv = V, Kv = -Av, (12) 

the existence and uniqueness of a local in time smooth solution [y, v) of 
the GNSB equations fllU|lip . for each smooth initial initial condition {y^,v^) 
given on the torus T'^, follow from standard theory on Euler and Navier- 
Stokes equations (for which we refer to [Li]). 
We say that {y, v) is a weak solution if: 

1) yit,x) and v(t,x) depends continuously on t with values in L'^{D,R'^) 
(with respect to the weak topology of L^); 

2) For all smooth time dependent vector fields w{t,x), z(t,x), with V -w = 0, 
we have: 

— y f ■ ti? = J[ev ■ {dt + {v ■ 'V))w — Kv ■ w + F{x, y) ■ w]dx, (13) 
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J y- zdx = J y-{dt + {v- V))z + G{x, y) ■ z)dx, (14) 



di 

3) The following energy inequality holds true: 

-— y ( e|t;p + lyp + J Kv ■ V dx < J [F{x,y) ■ v + G{x,y) ■ y] dx. (15) 

When K = —A, the existence of global weak solutions for the GNSB equa- 
tions follows from standard arguments a la Leray combined with the DiPerna- 
Lions theory on ODEs [Lil IDL] . They are unique in 2 space dimensions. In 
sharp contrast, as = / or i^' = 0, nothing can be said about global weak 
solutions. 

Concerning global smooth solutions, the existence theory is quite challenging, 
even in 2 space dimensions. Recently, Chae, Hou and Li |ChllHL] have proven 
that Navier-Stokes Boussinesq equations (IHllQl) (just called Boussinesq equa- 
tions in these papers) have global smooth solutions when d = 2 and K = —A. 
The same result can be readily extended to the GNSB equations ( ITOllTTl) . es- 
sentially because we assume the right-hand sides F{x, y) and G{x^ y) of each 
equation to be smooth functions of y and v with bounded derivatives up 
to order two. (Indeed, these assumptions are enough for a straghtforward 
adaptation of the proof of Theorem 1.1 in Chae's paper, through estimates 
(2.1 • • • 17) in [Chj . Some constants involved in these estimates have just to 
be modified to take into account the Lipschitz constants of F and G.) 
So, we can summarize all these results in the following Theorem, which is 
nothing but a straightforward adaptation of known results: 

Theorem 2.1 Assume the dissipative operator K to be of type [W\] . Then, 
the generalized Navier-Stokes Boussinesq equations / I70|f77]) admit, for any 
smooth initial condition, a unique local smooth solution. 
If K = —A, the GNSB equations admit at least a global weak solution {y,v) 
(in the sense ( f7gp^p5]) ) for any initial condition (|/°, t>°) in L"^. If d = 2, 
these weak solutions are unique. Furthermore, still for d = 2, the solutions 
are globally smooth for smooth initial conditions. 



2.2 Zero inertia limit of the GNSB equations 

By zero inertia limit of the GNSB, we mean the formal limit obtained by 
dropping the scaling factor e in front of the inertia terms in fllOlllip . Namely, 
in Eulerian coordinates, 

dty+{v-V)y = G{x,y), (16) 
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Kv + Vp = F{x,y), V-v = 0. 

We are able to make a rigorous derivation of the zero inertia limit when K 
is strictly dissipative {K = being so excluded): 

Theorem 2.2 Assume the dissipation operator K to be coercive in , namely 
K > a for some constant a > 0. Then the zero inertia equations [T^) admit, 
for any smooth initial condition, a local smooth solution, which is global if 
d = 2 and K = —A. This solution can be obtained as the limit, as e goes 
to zero, of the weak solutions of the GNSB equations ( flOlfTI]) . with the same 
initial condition. 

For the convergence, we use a simple energy method. Namely, given a 
weak solution {y' , v') to the GNSB equations ( |T3|14II15I) and a solution (y, v) 
of the HF equations, with same initial conditions (?/'^,t'°), we introduce 

e(t)= / {e\v'\' + \y-yr)dx (17) 

and try to get an estimate of form: 

|(e(t) + 0(e)) + \ j^^ K{v - v') ■ {v - v') dx < (e(t) + 0{e))c, (18) 

where c depends only on the limit {y,v), for any fixed finite time interval 
[0, T] on which {y, v) is smooth. From this estimate (fTSll . we immediately get 
that y — y' and v — v' are of order 0{^/e) in, respectively, L~([0, T], L2(T^)) 
and L^([0, T], L^(T°')), using the coercivity of K {K > a for some a > 0). 
So, we are left with proving ( JTSll . Notice first that, from equations ( JT5l) and 
f|T6|) . the following energy balances hold true: 

~Ji ^W\^ + b? )dx + J Kv' ■ v' dx < J [F{x, y') ■ v' + G{x, y') ■ y'] dx 

2^1 1^1^ '^^ ~'~ J ' ^ dx = J [F{x,y) ■ V + G{x,y) ■ y] dx. 

Since {y', v') and (y, v) are respectively supposed to be a weak solution of 
the GNSB equations and a smooth solution of the zero inertia limit f|T6|) . we 
also get 

-jtl yy''^'^ 
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[{{v ■ V)y) ■ y' - ((t/ ■ V)y) ■ y']dx - J [G{x, y) ■ y' + G{x, y') ■ y]dx 

{{v V)y - {v' ■ V)y) ■ {y' ~ y) dx 

+ j y) - y')) ■ {y - y') - G{x, y') ■ y' - y) ■ y]dx 

[using that / {w ■ V)y) ■ y) dx = for both w = v and w = v') 

< ci / (\v -v'\\y - y'\ + \y- y'\^) dx - [G{x,y') ■ y + G{x,y) ■ y]dx 



where ci depends on the Lipschitz constants of G (as a function of y and v) 
and V (as a function of x). Thus, adding up these three equahties, we get by 
definition ([T7|) : 

d r 

— e(t) + J {Kv' -v' + Kw) dx 



< Ci J {\v — v'\\y — y'\ + \y — y'p) dx + j [F{x, y') ■ v' + -F(x, y) ■ v] dx. 
From ([T3|) and f|T6|) . we also get: 

j F{x, y) ■ v' dx = j {Kv + Vp) ■ v' dx = J Kv ■ v' dx 



and 

d 



J F{x, y') ■ V dx = e— J v' ■ v dx — e J v' ■ {dt + v' ■ V)f dx + J Kv' ■ v dx 

= + r^it)) + e2(t) + I Kv'-v dx, 

where 



ri(t) = \/e y t>' ■ t> dx, 
'f2{t) = — a/c j v' ■ dtv dx 



and 

Notice that rf, r| and |e2| are bounded by C2e{t) (by definition (fTTl) ). where 
C2 depends on the Lipschitz constant of v as a function of both t and x. So, 
we have obtained: 

^(e(t) - v^ri(t)) + 1 ir(i;' - v) ■ {v' - v) dx 
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< ci y {\v - v'\\y - y'\ + \y - y'H dx + y/er2(t) + C2e{t). 

Using the coercivity of K {K > a > 0) and definition flTTl) . we find C3 
depending on the Lipschitz constants of F, G and v such that: 

ci J ~ "^1 b "~ y'\ + \y ~ y'f) dx < C3e(t) + 2 _/ -^('^' ~ ■ i^' ~ dx. 

This leads to: 

^(e(t) - v^ri(t)) + i 1 - v) ■ (v' - v) dx < (cs + c^)e{t) + v^r2(t), 

which leads to a differential inequality of the desired type, namely (ITSi) . since 
r\ + rl < 2c2e. Thus, the proof of Theorem 12.21 is now achieved. 



2.3 A Mechanical interpretation of the GNSB equa- 
tions 

In this subsection, we provide a mechanical interpretation of the GNSB equa- 
tions and their zero inertia limit. For this purpose, it is worth considering the 
GNSB equations fllO|lll) in "Lagrangian coordinates". Assuming the vector 
field V to be smooth enough, denoting by a G -D the position of a fluid parcel 
at time t = 0, we can recover its position a) G -D at later time t > by 
solving the ODE 

dtX{t,a) = v{t,X{t,a)), X(0,a)=a, WaeD. (19) 

Notice that, for each t, a E D ^ X{t,a) G -D is a measure preserving map 
as a consequence of the fact that f is a smooth divergence-free vector field 
parallel to dD. Let us also introduce: 

Y{t,a) = y{t,X{t,a)) e R"". (20) 

Then, the GNSB equations fllOflll) read in Lagrangian coordinates: 

e duX{t, a) + {Kv){t, X{t, a)) + (Vp)(t, X{t, a)) = F{X{t, a), Y{t, a)), (21) 

dtY{t,a) = G{X{t,a),Y{t,a)), 

where a E D ^ X{t,a) E D is Lebesgue measure preserving. Let us now 
provide a possible mechanical interpretation. We model the atmosphere (or 
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the ocean) as a continuous distribution of infinitesimal rigid balloons floating 
inside each of them having position X(t, a) at time t, with X(0, a) = a, 
and being attached with probability A (a) > to an an anchor with an elastic 
cable. Of course, to be a realistic model with real balloons, A should be a 
discrete probability distribution concentrated on a finite collection of points 
and the corresponding balloons should have a finite extension! Let us rather 
assume, for mathematical simplicity, that the balloons are just points and 
that A is a smooth nonnegative density function on D with unit mass. The 
cable corresponding to the balloon labelled by a is modelled by a (possibly 
non Hookean) spring with restoring force —k{C,, a) = k{—C,, a) G where 
^ G -R*^ is the elongation of the spring. Notice that k may depend on a. The 
location of the anchor attached to the balloon labelled by a is not necessarily 
fixed and denoted by Y{t,a) G W^. (We may also think of an aircraft, or 
a boat, or any kind of carrier instead of an anchor.) Notice that we do not 
require the anchor to be located in D. Neglecting any interaction between the 
fluid and both the anchors and the springs (which may not be very realistic), 
we obtain the following dynamical equation for each balloon 

edttX{t,a) + {Kv){t,X{t,a)) + (Vp)(t, X(t, a)) = (22) 

-X{a)k{X{t,a) - Y{t,a),a). 

(Observe that as A(a) = 0, the corresponding carrier Y(t, a) is just fictitious!) 
Let us consider the special case when the speed of each anchor is constant 
and given by: 

dtYit,a) = W{a) e R\ (23) 
Implicitly define a field y = {y, y) = y{t,x) & x D by setting 

y{t, X{t, a)) = Y{t, a), y{t, X{t, a)) = a 

(remember that a X{t,a) is supposed to be a diffeomorphism) . Noticing 
that 

{{dt + V ■ V)y)(t, Y{t, a)) = dt[y{t, X{t, a))] = {dtY{t, a), 0) 

= {W{a),0) = {W{y{t,X{t,a))),0) 

and going back to Eulerian coordinates, we recover the GNSB equations 
(llOflll) in the particular case: 



F{x,y) = -Xiy)ki^-y,y), G(x, y) = (iy(y), 0). (24) 
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(Notice that assuming A and the restoring force k to have bounded derivatives 
up to order 2 is not very reahstic! These assumptions are clearly made for 
mathematical convenience.) We may also consider the following variant of 
this mechanical model. Instead of prescribing their velocity by fl23|) . we may 
assume that the carriers are driven by a friction-dominated retroaction of 
type: 

r]dttY{t, a) + dtY{t, a) = —fi{a)k{Y{t, a) — X{t, a), a), 

where /i is a given nonnegative function. Dropping the inertia term [i] = 0) 
leads to the following law: 

dtY{t, a) = -fi{a)k{Y{t, a) - X{t, a), a), (25) 

In that case, we get (keeping unchanged the definitions of Y and y) again 
the GNSB equations ffTOllTTD with F and G given by: 

Fix,y) = -X{y)k{x - y,y), G{x,y) = {-ij,{y)k{y - x,y),0). (26) 

Let us finally consider a second variant where the Coriolis force is added to 
the model (rotating ocean or atmosphere). Neglecting the vertical extension, 
so that d = 2, and assuming the rotation vector to be perpendicular to the 
ocean and of unit length, the Coriolis force Jv = {—V2,vi) is completely 
absorbed by the pressure term and the fluid parcels are not sensitive to it. 
(Indeed, v = (fi,f2) being divergence free, Jv is a gradient and can be 
removed from the dynamical equation.) However we may think that the 
carriers are still sensitive to the Coriolis force. Thus, we get for them 

■r]dttY{t, a) + JdtY{t, a) = —^{a)k{Y{t, a) — X(t, a), a), 

instead of 0251) . neglecting a possible friction term. Neglecting the inertia 
term {rj = 0) leads to the balance equation: 

dtY{t, a) = Jfx{a)k{Y{t, a) - X{t, a), a) (27) 

(using that J"^ = —J). So we end up with a third version of the GNSB 
equations fllQ|lip . for which: 

F{x, y) = -\{y)k{x -y,y), G{x, y) = {Jf^miv - x,y),0). (28) 

To summarize this subsection, let us just say that the GNSB equations 
fllOflll) provide a rather flexible framework to describe the interaction be- 
tween an incompressible fluid confined in a rf— dimensional domain D (each 
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fluid parcel being labelled by its initial position a E D having position X(t, a) 
at time t), and a set of particles (also labelled by a and of position Y{t,a)) 
moving in the ambient space R^. The unusual feature of the resulting mod- 
els is that the interaction is pairwise {X{t,a) interacts only with Y{t,a) for 
the same label a), except for the mediation by the pressure field p{t, x) wich 
preserves a — > X{t, a) in the class of Lebesgue measure preserving maps of 
D at each time t. 

3 The Generalized Hydrostatic Boussinesq equations 

In this final section, we investigate the most degenerate version of the GNSB 
equations fllOflll) . where we neglect not only the inertia terms but also the 
dissipative operator K. Thus we are left with the strange looking system: 

F{x,y) = Vp, V-v = 0, (29) 

dty + {vV)y = G{x,y), (30) 

that we call Generalized Hydrostatic Boussinesq (GHB) equations (by seing 
fl29|) as a generalization of the hydrostatic balance in Convection Theory). 
Let us concentrate on the simpler case when m = d and 

Fix,y) = y-x (31) 

(which corresponds to cables modelled by Hookean springs, according to the 
mechanical interpretation of subsection 12.31) . Thus, (!29l) just reads: 

y = x + Vp{t,x). (32) 

In Lagrangian coordinates, the Generalized Hydrostatic Boussinesq equations 
(15011521) become: 

Y{t, a) = X{t, a) + (Vp)(t, X{t, a)), (33) 

dtY{t, a) = G{X{t, a), Y{t, a)), (34) 
where, for alH, a G -D ^ X{t, a) G -D is a measure preserving map. 

3.1 Formal derivation of some optimal transport mod- 
els from the GHB equations 

We claim that several models involving optimal transport and the Monge- 
Ampere equation correspond to these GHB equations. In particular, we 
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consider the following generalization of the "semigeostrophic (SG) equations" 
[miCNPllBBllCGPl lL^: 



dtp + V- {pw) = 0, (35) 

w{t,x) = BV!f, (36) 

det{I + DV) = P, (37) 

where B is a d x d constant matrix and D'^if is the "Hessian" matrix, made 
of all second order derivatives of fit,x) with respect to x. This system, that 
we call Generalized Semi-Geostrophic (GSG) equations involves the "fully 
nonlinear" Monge-Ampere (MA) equation (!37l) which, requires, in order to 
be of elliptic type, the convexity condition: 

I + D^ip{t,x) > 0, (38) 

in the sense of symmetric matrices, for each time t. The 2D SG equations 
|Ho[ IGNPl IBB[ IGGP[ ILoj just correspond to the special case when d = 2 and 
B is the rotation matrix of angle 7r/2. The case when B is just a number can 
be related to drift-diffusion and Keller-Segel type models |CMPS] for which 
the MA equation is replaced by the linear Poisson equation 

A<^ = p-1, (39) 



which can be seen as a linear approximation of the MA equation fl37j) as 
if is small. The drift-diffusion case corresponds to fl35ll36|39p with B > 0. 
The simplified version of the Keller-Segel model |KSj treated by Jager and 
Luckhaus [JL] corresponds to S < (with an additional diffusion term for p 
in equation (l35l)). 



Let us now show that a solution of the GHB equations (133II34I) corresponds 
to a solution of the GSG equations fl35ll36ll37|l38!) . in a suitable sense. For 
this purpose, in order to use the Polar Factorization Theorem II. 3[ we make 
the following a priori assumptions for each time t: 
Al: The map Y{t,-) is non degenerate, 

A2: The map x E D ^ x + Vp(t, x) G has a convex potential. 

These assumptions mean that fl33|) defines the polar factorization ofY{t, ■) 
where X{t, ■) is measure preserving and x E D ^ x + Vp{t, x) has a convex 
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potential. According to the Polar Factorization Theorem 11.31 the measure 
preserving factor X{t, ■) can be written: 

X(t,a) = (V$)(t,r(t,a)), 

where x) is convex and Lipschitz continuous in a; G W^, or, equivalently, 

X{t, a) = Y{t, a) + (V<^)(t, Y{t, a)), (40) 

where f(t,x) = x) — |xp/2. Since Y is supposed to be non degenerate (by 
assumption Al), there is a nonnegative Lebesgue integrable "density field" 
p{t, x) such that 

/ f{x)p{t,x)dx= f f{Y{t,a))da, (41) 

for all suitable functions /. Thus: 

/ f(x + V(p(t,x))p(t,x)dx = [ f(Y(t,a) + (Vip)(t,Y(t,a)))da 

(by definition (jUj) of p) 

= I f{X{t,a))da= / f{x)dx 
Jd Jd 

(thanks to ( HOj) and because X(t, ■) is Lebesgue measure preserving). So, we 
have obtained 

/ f(x + V(p(t,x))p(t,x)dx= I f(x)dx, (42) 
jR'i Jd 

for all compactly supported continuous function /. This, combined with 
the assumption that x + Vp(t, x) has a convex potential can be seen as a 
weak form of the Monge- Ampere equation fl371) combined with the ellipticity 
condition (I55|) . Next, using fHUl) . we can write as 

dtY{t,a)=w{t,Y{t,a)), 

where w is the vector field defined by: 

w{t,x) =G{x + V(p{t,x),x), (43) 
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which is nothing but a generahzation of equation (]36ll . Next, we get for all 
smooth compactly supported function / on R^: 

^ f{Y{t, a))da = jjyf){Y{t, a)) ■ w{t, Y{t, a))da, 



which means, in terms of p defined by (j?T 
d 



I f{x)p{t, x)dx = / V f{x) ■ w{t, x)p{t, x)dx, 
dt J jR'i 

that is (135|) in a weak sense. So, we have fully recovered the GSG system 
fl35f36ll37|l38|) . in a suitable weak form, from the GHB equations ( !33|MI) under 
Assumption Al,A2. In addition, equation (136|) can be replaced by the more 
general relation 



3.2 A global existence theorem for the GHB equations 

We are now going to introduce a suitable concept of solutions of fl33|34p . by 
assuming a priori that, in (1551) . for each time t, the map x G -D — > x+Vpit, x) 
has a convex potential and, therefore, is the unique rearrangement Y*{t,-) 
of Y{t, ■) in the class C of map with a convex potential, as in Theorem 11.21 
(This corresponds to assumption A2 in the previous subsection.) Therefore, 
we can just write (j33ll as: 

Y{t,a) = Y*{t,X{t,a)). 
From ( 134|) . we get (at least formally) that: 

^ f{Y{t, a))da = J^iVf){Y{t, a)) ■ G(X(t, a), Y{t, a))da. 

for all compactly supported function /. Thus, 

jJ^f{Y\t,X{t,a)))da = 

= I (Vf)(Y*(t,X(t,a))) ■G(X(t,a),Y*(t,X(t,a)))da. 

JD 

Now, we can factor out X{t, a) and get a set of self-consistent equations for 
Y*, namely: 

^ / f{Y*{t,a))da= f {Vf){Y*{t,a))-G{a,Y*{t,a))da, 
dt Jd Jd 
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without loss of information for Y*. This suggests the following concept of 
solution for the GHB equations ( 133)1341) : 



Definition 3.1 Assume that G is continuous and satisfies 

sup < oo. (44) 

We say that Y* E C°{[0,T], L'^{D, R'^)) is the "convex rearrangement" (CR) 
solution to the GHB equations ^3^3^, if: 

1) F*(t, ■) belongs to the set G of all maps with convex potenial, for all 

te[o,T], 

2) For all compactly supported G^ function f on R^, we have: 

/ /(F*(t, a))da = I {Vf){Y*{t, a)) ■ G(a, F*(t, a))da. (45) 

dt JD JD 

This concept yields the following global existence theorem (without unique- 
ness) for all initial conditions in L^: 

Theorem 3.2 For each initial condition Y^ G L^{D,R'^), there is at least 
one CR-solution Y*{t,a), in the sense of Definition \3.1\ 
such that Y*{t = 0, ■) = {Y^ . 

This solution can he obtained as the limit in C°(L^) as h of a time 
discrete approximation Y^{t,a) defined, first at discrete times t = nh, by: 

Y''{nh + h,a) = [Y''{nh,a) + hG{a,Y''{nh,a))]*, n = 0,1,2,- ■■ (46) 

(where * denotes the rearrangement operator as in Theorem \l.S\} and, then, 
linearly interpolated in t. 

Proof 

To get the existence result, it is enough to show the convergence of the time 
discrete approximation Y^. First, we observe that, Y^{t, ■) is valued in G, 
(the class of maps with convex potential) for all time t. (This is true by 
definition for discrete times t = nh and preserved by linear interpolation 
since C is a convex cone.) Next, we deduce from fHBl) and assumption 



I j^\Yh{nh + h,a)\'^da<hc+{l + hc)J j^\Yh{nh,a)\'^da, (47) 
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for some constant c depending only on G and D. We also get, for all com- 
pactly supported function /: 

[f{Y^\nh + h, a)) - f{Y^\nh, a))]da = (48) 



= h [ C{Vf){Y^{nh, a) + heG{a, Y^{nh, a))) ■ G(a, Y^{nh, a))d9da, 

JD Jo 

which can be bounded by: 

he sup ^-^^^ J^{l + \Y\nh,a)\^)da, (49) 
where c depends only on D and G. 

So, from (H7j) . we first see that Y^{t, ■) is bounded in uniformly in t e [0, T] 
and h g]0,T] by some constant R. Therefore Y^ is uniformly valued in Cr 
the set of maps with convex potential and norm bounded by R, which is 
a compact subset of L^. 

Next, we deduce from ( H8f49l) that, for each fixed function / such that 



|V/(x) 



\x\ 



< +CX), t [ f{Y^{t,a))da 

JD 



is Lipschitz continuous on [0,T], uniformly in h. Since Theorem 11.21 asserts 
the continuity of the map /i — > y*, we get that Y^{t, ■) is uniformly equicon- 
tinuous from [0,T] to L^. Then, we deduce from the Ascoli-Arzela theorem 
that the set of all time-discrete approximations Y^{t,a), for < /i < T, is 
relatively compact in C°(L^). Thus, there is a sequence of time steps h for 
which Y^ converges to some limit Y* in Gf{Ll), which is necessarily valued 
in Gr, and therefore in G. We also easily get fHS]) by letting h go to zero in 
So, the proof of Theorem 13.21 is now complete. 



Remark: continuous dependence as d = 1 

In the very special case d = 1, the rearrangement operator * is well known 
to be non expansive in L^: 

/ \y*{x) — z*{x)\'^dx < / \y{x) — z{x)\'^dx, 

JD JD 
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for pair (y, z) of applicatfons from D to R. It foUows that two CR solutions 
y* and z*, obtained as limits of the time-discrete approximations ( H6ll . with 
respective initial condition and z^, must satisfy: 



where c depends only on D and G. 

4 Optimal transport and Magnetic Relaxation 

In this section, we discuss a natural "stringy" generalization of the AHT 
equations (Pl^ (discussed in the first part of the paper) and estabhsh a hnk 
with the Arnold-Moffat model of Magnetic Relaxation (see |AK[ [Mol IMo2t 



4.1 The AHT model as a gradient flow 

Using Lagrangian coordinates, we deduce from (J3f^ . in the case when the 
dissipation operator is = 1: 



where X{t, ■) belongs to the set MPM{D) of all (Borel) Lebesgue measure 
preserving maps of D. (These equations can be either derived directly from 
(13TO or obtained from the GNSB equations written in Lagrangian coordi- 
nates ( 12T]) . by setting e = 0, K = I, F{x,y) = y — x and G{x,y) = 0.) 
As explained in [AHT], in slightly different words, the AHT equation (!50i) 
formally corresponds to the gradient flow of the "energy" 



on the "manifold" of all X G MPM{D) for the metrics. Let us recall 
that, as seen in the first section, minimizing this energy is equivalent to solve 
an Optimal Transport problem. The gradient flow structure can be easily 
understood by considering the standard time discretization of such a gradient 
flow. Let h > he a time step and let us denote by X^{t,a) the discrete 
approximation of X{t, a) at discrete time t = nh, n = 0, 1, 2, 3, ■ ■ ■. At t = 0, 
we set X^{0, a) = a and, for t = nh, n = 1, 2, 3, ■ ■ ■, we require X^{t, ■) to be 
a minimizer among all X G MPM{D) of the following functional: 




l5H IVMn iNi]l 



dtX{t,a) + {Vp){t,X{t,a)) = y^{a)-X{t,a) 



(50) 




(51) 




(52) 
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or, equivalently, 

/ |X(a)-^'("'-''"^+^°("^'prfa (53) 
Jd^ ^ ^ l + h ' ^ ^ 

(after rearranging the squares). 

Thus, assuming a priori that X^{nh — h, ■)+y'^h is non degenerate and using 
Theorem 11.31 this exactly means: 

X\nh -h,-) + y% = y*o X\nh, ■) 

where y* is a map with a convex potential. If we write this map as 

X ^ X + h{Vp^){nh, x), 

we get 

X\nh -h,-)+ y% = X^{nh, ■) + h{Vp^){nh, X^{nh, ■)), 

which can be seen just as a finite difference approximation of equation (ISOll 
as h 0. This is enough to interpret, at least formally, equation (1501) as the 
gradient flow of energy ( ISTl) on the "manifold" MPM{D). 

4.2 A stringy generalization of the AHT model 

This analysis suggests the possibility of more complex models based on 
similar ideas. A rather natural idea amounts to consider, instead of the 
"manifold" MPM{D), the "manifold" of "strings" valued in MPM{D): 
s e [0, 1] X(-, s) e MPM{D), with fixed end values, say 

X{t,a,s = 0) = X~{t,a), X(t, a, s = 1) = X+(t, a). (54) 

Then, we may think of the gradient flow of the following string energy: 

- /V \dsX{a,s)\'^dads. (55) 
2 JO Jd 

We claim that the resulting equation read: 

dtX{t,a,s) = (9f,X(t,a,s) + (Vp)(t, X(t, a, s), s), (56) 

where X{t,-,s) is valued in MPM[D) and end point conditions fl5^ are 
enforced. To get this system, as we did before for the AHT model, we define 
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a time discrete approximation X^{t,a, s), by setting X^{0,a,s) = a and 
asking, for t = nh, n = 1,2,3, ■■ ■, X^(t, -, ■) to be a minimizer among all 
curves s e [0, 1] X{-, s) e MPM{D) of the functional: 

Jd Jo 2h Jd Jo 2 

The formal optimality condition reads: 

X'^{nh, a, s) — X^{nh — h, a, s) 

h " 

= dlX\nh,a,s) + {Vp''){nh,X''{nh,a,s),s), 

for some scalar function p^. So, we formally obtain, as /i — 0, the desired 
equation fl56p . Equation (1561) has an interesting interpretation, obtained 
by assuming a priori that a & D ^ X{t,a,s) is a smooth orientation and 
measuring preserving diffeomorphism of D for each {t, s). Then we introduce, 
for each (t, s), two divergence free vector fields parallel to the boundary dD, 
namely v(t, x, s) G R'^ and b{t, x, s) G R^, defined by: 

f (t, a, s), s) = a, s), b{t, X{t, a, s), s) = dsX{t, a, s). (58) 

Then, we get from (l56il : 

v = dsb+{b- V)6 + Vp, (59) 

while, from fl58|) . we get the compatibility condition 

dtb + {v ■ V)6 = (b ■ \/)v + dsv (60) 

(by writting d^^X = d%X), to be added to the divergence free constraints 

V-i; = V-6 = 0, v//dD, b//dD, (61) 

and the boundary conditions at s = and s = 1 induced by (I54p . namely: 

v{t,x,s = 0) = v~(t,x), v(t,x, s = 1) = v^(t,x), (62) 

where and v~ are prescribed. When the fields v and b do not depend on s, 
we get the Magnetic Relaxation model discussed by Moffatt in |Moj (see also 
[AKl IMo2l [gel IVMll [Nij). As t -> +oo, we expect, at least for a large class 
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of initial conditions, the solution of equations fl59ll60|6ip to converge toward 
an equilibrium, for which v = and b = b{x,s), p = p{x,s) are solutions 
to the Euler equations |AK[ IMPj [s acting as the time variable and b as the 
velocity field): 

dsb + {b ■ V)6 + Vp = 0, (63) 

V ■ 6 = 0, b//dD. (64) 

Of course, we are far from being able to provide any rigorous proof of this 
conjecture. 



4.3 The "cross-Burgers" equation 

In the case when D is the unit ball. The Magnetic Relaxation equations 
(I59|60ll6ip admit special solutions (6, f , Vp) which are linear in x: 

b{t,x, s) = B(t, s)x, v(t,x, s) = V(t, s)x, Vp(t,x, s) = G(t, s)x, (65) 

where B, V are skew-symmetric matrices, while G is a symmetric matrix, all 
depending only on (t,s). (Notice that the fields b and v are automatically 
parallel to the boundary dD since D is the unit ball.) The resulting equations 
for B, V and G are: 

V = dsB + B'^ + G, (66) 
dtB + [V, B] = dsV. (67) 
Since B^ is a symmetric matrix, equation ( l66i) reduces to: 

V = d,B. 

Thus, we get a single equation for B: 

dtB + [dsB, B] = dlB, (68) 

where [A, B] denotes the skew product AB — BA. In the special case d = 3, 
B can be identified as a 3- vector and [■, ■] as the cross product x in i?'^, which 
leads to: 

dtB + d,Bx B = dlB. (69) 

that we could call the "cross-Burgers" equation. This equation admits inter- 
esting special solutions, such as: 

B{t, s) = (ait) cos s , a{t) sin s , /3(t) — 1) 
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where a > and (3 are solutions to: 




or, equivalently, 



+ exp(2A) = 0, 



df" 



where A = log(Q;). 
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